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ABSTRACT 

We present a unified picture for the evolution of star clusters on the two- 
body relaxation timescale. We use direct N-body simulations of star clusters in 
a galactic tidal field starting from different multi-mass King models, up to 10% 
of primordial binaries and up to N tot = 65536 particles. An additional run also 
includes a central Intermediate Mass Black Hole. We find that for the broad range 
of initial conditions we have studied the stellar mass function of these systems 
presents a universal evolution which depends only on the fractional mass loss. 
The structure of the system, as measured by the core to half mass radius ratio, 
also evolves toward a universal state, which is set by the efficiency of heating on 
the visible population of stars induced by dynamical interactions in the core of 
the system. Interactions with dark remnants (white dwarfs, neutron stars and 
stellar mass black holes) are dominant over the heating induced by a moderate 
population of primordial binaries (3-5%), especially under the assumption that 
most of the neutron stars and black holes are retained in the system. All our 
models without primordial binaries undergo a deep gravothermal collapse in the 
radial mass profile. However their projected light distribution can be well fitted 
by medium concentration King models (with parameter Wq ~ 8), even though 
there tends to be an excess over the best fit for the innermost points of the surface 
brightness. This excess is consistent with a shallow cusp in the surface brightness 
(fj, ~ R~ u with v ~ 0.4 — 0.7), like it has been observed for many globular clusters 
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from high-resolution HST imaging. Generally fitting a King profile to derive the 
structural parameters yields to larger fluctuations in the core size than defining 
the core as the radius where the surface brightness is one half of its central value. 
Classification of core-collapsed globular clusters based on their surface brightness 
profile may thus fail in systems that appear to have already bounced back to lower 
concentrations, particularly if the angular resolution of the observations is limited 
and the core is not well resolved. 

Subject headings: Stars: luminosity function, mass function - Galaxy: globular 
clusters: general - Methods: n-body simulations - stellar dynamics 



Introduction 



Until a few years ago, the standard picture of globular clusters stellar population and 
dynamical evolution emerging from theoretical and observational studies was yielding a con- 
sistent and well understood framework. Globular clusters were thought to be 'simple stellar 
population', composed of stars with the same age and chemical composition. Globular cluster 
dynamical evolution had been thoroughly investigated by a large number of numerical stud- 
ies showing that clusters surviving the early expansion triggered by primordial gas expulsion 
and mass loss due to stellar evolution evolve, due mainly to two-body relaxation, toward 
core collapse and higher central concentration s while losing stars th rough the boundary set 
by the tidal field of their host galaxy (see e.g. iHeggie fc Hut II2003I ). 



In the last few years, however, a wealth of observational data mostly from exquisite 
space based observations have revealed the actual complexity of globular cluster dynamics 
and stellar populations. It is now clear that there is a close link and interplay between 
dynamical evolution and the stellar content of clusters, the structure of clusters and the 
abundances of exotic objects (such as blue stragglers, X -ray sources, pulsars etc.) (see e.g. 
Belczvnski etallbood : Ishara & Hurlevlbood : iHufaod ). 



Recently, a number of photometric and spectroscopic observations have also shown evi- 
dence of the presence of multiple stellar populations and challenged the commonly held view 
according to which globular cluster are 'simple stellar populations'. It appears that essen- 
tially any globular cluster that has high quality photomet ric and spectroscopic data avail- 
able presents evidence against a single star-formatio n burst (IGratton et al.ll2004j : iBedin et al. 
2004 : iPiotto et al.ll2005l 120071 : ICarretta et al.l 120091 ) . The presence of multiple stellar popu- 
lations h as major implication s for the formation and the early evolution of globular clusters 



^see e.g. 



D'Ercole et all 120081 ). 
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Furthermore, our understanding of the late-ti me evolution of globular clusters has been 
called into question. A recent observational study iDe Marchi et al.l (120071 ). focused on the 
stellar mass function of Galactic globular clusters, found that clusters with lower concen- 
tration index of the surface brightness profile tend to have flatter stellar mass functions. 
The observed trend is at odds with what is expected from the standard dynamical evolution 
scenario according to which as clusters evolve toward core collapse and higher concentra- 
tions they lose mass and the preferential loss of low-mass stars flattens the stellar mass 
function; dynamically older clusters shou l d be m ore centrally concentrated and have flatter 
stellar mass functions. iBaumgardt et al. I (120081 ) have suggested that strong primordial mass 
segregation might be required to explain the fl at stellar mass function of some of the low- 
concentration cluster in the observed sample of IDe Marchi et al.l (120071 ). However a strong 
primordial mass segre gation might lead to t he rapid dissolution of the cluster due to stellar 
evolution mass- losses (jVesperini et al.ll2009i ). 



High-res olution pho tometry of the center of globular clusters has also revealed that 
the standard iKingl (119661 ) models do not capture the details of the surface b rightness profile, 
which often shows the presence of shallow cusps within a relatively large core (INoyola fc Gebhardt 
20061 ). The shallow cusps contrast with the classic expectation that a strong, isothermal cusp 
develops as a result of core-colla pse. The presence of a central IMBH is a p ossible explanation 
for the presence of those cusps (IBaumgardt et al.ll2004j ; iTrenti et al.ll2007l ). but this interpre- 
tation a ppears in contra st with the lack of strong evidence for IMBHs in globular clusters 



l^e.g. see lGill et al.ll2008l . and references therein). NGC2298, for exampl e, is one of the clus- 
ters sh owing a central excess in the surface brightness profile from t he iNovola &: Gebhardt 



( 120061 ) data, but a central IMBH is ruled out at high confidence level ( jPasquato et al. 



2009 



The goal of this paper is to study in detail the relation between a cluster structural 
properties and the dynamical phases of its evolutionary path as well as the relation between 
the properties of a cluster stellar mass function and its structure. We focus our attention 
on the long-term evolution of star clusters driven by the effects of two-body relaxation. We 
analyze our numerical models following the procedures used in the observational studies 
of the structure and stellar content of star clusters and by directly comparing our results 
with observational data we aim at clarifying their interpretation in the context of dynamic 
evolution of relaxed stellar systems. 

We have carried out a survey of multimass simulations of star clusters evolving in a tidal 
field and starting from a broad range of different initial conditions spanning different initial 
density profiles, primordial binary fractions and initial stellar mass function. This paper 
is organized as follows. In Section [2] we present our sample of simulations. In Section [3] 
we introduce our framework for analyzing simulations consistently with the observational 
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derivation of structural quantities. In Section H] we presents the results of our analysis, 
concluding in Section O 



Numerical framework 



2.1. Past investigations 

Numerical investigations of the dynamics of globular clusters are computationally very 
challenging because of the need to resolve a large dynamic range in timescales, from a few 
hours typical of a tight binary to the billion of years over which the system evolves due to 
two-body relaxation. In addition, the computational complexity of a run carried out for a 
constant number of relaxation times scales as i V 3 / log(iV), b ecause the relaxation time t r h 



increases with the number of particles N (e.g see ISpitzerlll987l ) . The computational resources 



required are further increased by the presenc e of a significant number of binaries, which can 



be up to 50% of the mass of the core (e.g., see iRubenstein fc Bailynlll9971 ; lAlbrow et al.ll2001 



Pulone et al.l 120031). al though some clusters, like NGC 6397, might also be almost binary- 
free ( iDavis et al.ll2008r). Hard binar i es have an important dynamic effect on the evolution of 



the c ore size flMcMillan et al.lll990l . 1 199 it IVesperini & Chernofllll994l ; iHeggie. Trenti & Hut 



20061 ). In the past, numerical simulations have b een performed either using approximate 



algorithms such as Monte Carlo methods (e.g., see iGao et al.lll99ll ; iFregeau fc Rasiol 12 007) 



or dir ect N-body simulations with a modest number of particles (N ~ 10 3 in iMcMillan et al. 
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19901 and IHeggie Aarsethl Il992l) until a recent improvement of one o rder of magnitude 



Heggie. Trenti fc Hutl 120061 : iTrenti. Heggie fc Hutl 120071 : iTrenti et all 12007 



20081 ; lHurleyil2007l ). These investigations highlighted two key results: (i) the presence of a 
significant population of primordial binaries drives the evolution of the core radius toward a 
value that is 5 — 10% of the half mass radius and (ii) the ratio of the core to half mass radius 
i r c/ r h) does not depend either on its initial value or on the primordial binary fraction /, 
once a few relaxation times have passed and provided that / > 0.1. In addition, runs with 
a galactic tidal field show the central concentration parameter c decreases toward the end 
of the simulation, due to the combined effect of binaries , that keep a large core, and of the 
evaporation, that progre ssively reduces the tidal rad ius (ITrenti. Heggie &: Hutll2007l ). How- 
ever, the simulations of ITrenti. Heggie fc Hutl (120071 ) were obtained using simplified initial 
conditions with equal mass particles and therefore cannot be directly applied to the study 
of the mass function evolution or analyzed consistently with observations. In this paper we 
extend our earlier work to include a realistic mass spectrum, as discussed below. 
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2.2. Our simulations 



T he simulation fram ewo rk that we use here has been extensively described in lTrenti. Heggie fc Hut 



( 120071 ) , I Gill et al.l (120081 ) and iPasquato et al.l (120091 ) . In short, we fo llow the evolution of star 
clusters using the direct summation code NBODY6 (jAarsethll2003l ). which guarantees an ex- 
act treatment of the multiple interactions between stars by employing special regularization 
techniques, without resorting to the introduction of softening. 



Our N-body simulations are carried out in natural (dimensionless) units ( iHeggie fc Mathieu 

19861 ). in which: 

G = M T = —AEt = 1, (1) 

where is the total mass of the system and Ej- the total energy (potential plus kinetic). 
The corresponding unit of time is: 



t d 



(GMj 



\5/2 



(4£ T ) 3 / 2 



(2) 



which approximately corresponds to the orbital period at the half-mass radius. In these 
units, the relaxation time can be written as: 



t 



rh 



0.138Nr 3 h /2 
log(0.11N)' 



(3) 



where ~ 1 is the radius at which the relaxation time is defi ned, that is the h alf-mass 
radius. Equation [3] can be translated in physical units as follows (IDjorgovski Ill993l ): 



8.9 x 10 5 yr 1M 

t r h — - ~ - - x ~ r X 



log (0.11 AT) (m,) \M e 



M 7 



1/2 



x 



lpc 



3/2 



(4) 



where (m,) is the average mass of a star in the system (including dark remnants such as 
neutron stars). 

The individual particle mass is drawn from an initial mass function appropriat e to study 
the la te evolutionary stage s of star c luster s, wh en stars are > 1 Gy r old. As in iGill et al. 



(120081 ) we consider either a ISalpeterl (119551 ) or a iMiller fc Scald (Il979l ) initial mass function, 
that is: 

£(m) cx m a , (5) 

with a = -2.35 and m e [0.2 : 1OO]M for the Salpeter IMF, while for the Miller & Scalo 
IMF the power-law slope is the following: a = —1.25 for m G [0.2 : 1]M , a = —2.0 for 
m G [1 : 2]M , a = -2.3 for m G [2 : 1O]M , and a = -3.3 for m G [10 : 1OO]M . 
Before starting the run, an instantaneous step of stellar evolution is taken to evolve the 
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mass function to a 0.8 M turnoff using the iHurley et al.l (120001 ) evolutionary tracks. In 
our standard models we have a 100% retention fraction of dark remnants. We also run two 
simulations with a 30% retention fraction of neutron stars and stellar mass black holes (see 
Table [T]). In all simulations no kick velocities are given to the remnants. Our idealized 
treatment of stellar evolution is based on the approximation that most of the relevant stellar 
evolution occurs on a timescale shorter than a relaxation time. It is thus appropriate to 
study the evolution of old globular clusters (age of ~ 10 Gyr). In fact, most of the impact of 
stellar evolution on the dynamics of a star cluster happens within the first few hundred Myr, 
when the most massive stars lose a significa nt fraction of mass and consequen tly contribute 



to a global expansion of the system (e.g. see iHurleyl 120071 ; iMackey et al.ll2008l ). Later in the 



life of a star cluster, two- body relaxation tends to erase the memory of the initial density 
profile and concentration (iTrenti. Heggie &: Hutl 120071 ; iTrenti et al.ll2008[ ). 



The initial distribution in the position- velocity phase space is that of a King model with 
scaled central potential Wq = 3, 5, 7. Our standard models have a primordial binary ratio 
of / = Nb/(N S + Nf,) between and 0.1, with N s and iVj, being the number of singles and 
binaries respectively. In addition, we als o consider a run with a cent ral Intermediate Mass 



Black Hole (discussed in more detail in iGill et al.l 120081 ). Following ITrenti. Heggie fc Hut 



(120071 ). our initial distribution of the binaries' binding energies is flat in log scale in the 
range from e m j„ to 133e m j n , with e m i n = (^(0))cr c (0) 2 . Here <r c (0) is the initial central 
velocity dispersion and (m(0)) is the average stellar mass at t = 0. For a typical velocity 
dispersion cr c (0) = 10 km/s, the semi-axes of binaries in the initial conditions are smaller 
than 10 AU. This choice is motivated by the fact that wider and therefore softer binaries, 
which are likely created in young star clusters, are destroyed within a fe w dynamical times 
and thus would have no effect on the late-time evolution of a star cluster (lHeggidll975l ). The 
stellar mass and type of binary members are drawn at random from our initial mass function 
and thus binaries containing either one or two dark remnants are possible. 

The models considered in this paper are tidally limited. The Galactic tidal field is 
treated as that of a point mass, and the tidal force acting on each particle is computed 
using a linear approximation of the field. The tidal cutoff radius of the simulation (which we 
indicate as {r t } s , where the curl-bracket with the suscript "s" is used to denote the definition 



of a simulation based quantity, see Section [3]) is defined as (ISpitzerlll987l ): 

M T 



{rt}] 



-R, 



Gah 



(6) 



SM Gal 

where Mc a i is the galaxy mass, and Real is the distance of the centre of the star cluster 
from the centre of the galaxy. We assume Real is constant, i.e. the cluster describes a 
circular orbit. The initial value of {r t } 3 is fixed to be self-consistent with the King profile 
used (for example, the Wq = 7 king model has {r t } s = 7.02 in N-body units). Particles are 
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rem oved from the system when they reach a distance 2 {r t } s from the center of the cluster 



see 



Trenti. Heggie &: Hutl 120071 for a full description of the tidal field treatment). 



We explore a large parameter space in the initial conditions varying the initial King 
model concentration, the Roche lobe filling factor, the initial mass function and the pri- 
mordial binary fraction, as summarized in Table [TJ The NBODY6 code has been run with 
its Graphic Processing Unit extension on the NCSA Lincoln Cluster. Each simulation was 
assigned to a computing node with 8 cores and two Tesla C1060 GPUs. We measured a 
sustained computational performance above 0.5 Tflops for a single node, that allowed us to 
simulate systems with up to iV = 65536 (or N = 64k in a compact notation where lk = 1024) 
particles until either t = 8000 was reached (corresponding to several initial relaxation times 
- t > 16 t rh (0)) or at least ~ 80% of the initial mass is lost due to evaporation of stars. 



3. Structural parameters of globular clusters: definitions in simulations and 

observations 

Both space-based photometry with HST and adaptive optics from the ground have 
greatly improved our knowledge of Galactic globular clusters. In fact, the angular resolution 
currently available is such that it is possible to resolve faint individual main sequence stars 
at the center of most Galactic globulars. One of the most striking examples is the globular 
clusters NGC 2298, where HST ph otometry reaches a completeness greater than 5 0% for 



0.2 M at the center of the system (IDe Marchi fc Pulond 120071 ; iPasquato et al.ll2009l ). 



A large fraction of globular clusters have been observed with HST (e.g. see lSarajedini et al 



20071 ). thus acquiring detailed information on their current stellar mass function and their 
central surface brightness profiles. Here we discuss the definition and determination of struc- 
tural parameters in both simulations and observations, with the goal to present a consistent 
analysis of our numerical investigation. Typically, structural properties in numerical simu- 
lations are based on a complete knowledge of the system with a tri-dimensional mass-based 
approach (e.g. the 3D half-mass radius of the system). Observations rely instead on pro- 
jected two-dimensional data, that are light based (e.g. the 2D half-light radius). In addition, 
observations are affected by finite angular resolution, limited signal-to-noise and field of view. 
To ensure self-consistency between quantities in simulations and observations, we developed 
a framework to "observe" a simulation snapshot, from which we obtain "observed" structural 
parameters. For reference we also provide structural parameters using the usual theoreti- 
cal/numerical definitions. Quantities that derive from the "observation" of a snapshot have 
no markings. 
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We adopt instead a curl-bracketing with an underscore "s" ({} s ) to denote the definition 
of a quantity based on current practice in numerical simulations (the "s" subscript stands 
for simulation). 



3.1. Observation of a simulation snapshot 



The positions and velocities of all particles in our simul ations are saved every 1 code 
time units (one time unit is about one dynamic time; see iHeggie fc Mathieul Il986l ). We 
then proceed to construct a synthetic observation by first ignoring all particles except main 
sequence stars. As our simulations only have an instantaneous step of stellar evolution, we 
don't consider stars along the giant branch. This choice is appropriate for comparison with 
star counts maps at high angular resolution, su ch as HST observation s, because giants are 
typically masked out of the analysis (e.g. see iDe Marchi et al.l 120071 ). Focusing on main 
sequence stars allows us to significantly reduce the surface brightness profile fluctuations 
that arise from the small number of luminous giants that would otherwise dominate the 
light profile. 

We then select a random direction and project the main-sequence stars positions creating 
a 2-dimensional map. W ithin this map we search for the center of the system using the 
Casertano &: Hutl (119851 ) density-center method. 



A circular surface brightness profile is then constructed using Int[y/ 'Nms] annuli (where 
N M s is the number of main sequence stars) each containing ~ \/N MS sources. From this 
profile we define the half-light radius in projection r^p as the radius containing half of 
the total luminosity. The King concentration parameter c = l o gm (r t/r r ) is determined 



by fitting the surface brightness profile with a single-mass iKingl (119661 ) model. We use a 
precomputed King profile table with uniform spacing in Wo (A Wo = 0.1). The fit is carried 
out imposing that the total luminosity and half-light radius of the fitting model match those 
of the simulation snapshot within a < 1% deviation. From Wq and thp we derive r t and r c . 
As an alternative definition for the observed core radius we consider the radius at which the 
surface brightness profile reaches half of its central value and we indicate this quantity with 



A'- 



To construct the global mass function of the system, we proceed to identify clo se pairs of 



stars w hose light would be b l ended together, similarly to the procedure described in lGill et al. 



(120081 ) and iPasquato et al.l (120091 ). Our main aim is to treat binaries in the simulations 
consistently with observations. We capture also line-of-sight superpositions and therefore 
crowding. For reference we set an HST-like resolution of 0.05 arcsec and we consider the 
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system, of assumed half-light radius 3 pc, at a distance of 5 Kpc from the Sun. This implies 
that we consider blended those stars separated by less than 250 AU, which includes all our 
primordial hard binaries. In addition, we get a small number of apparent binaries, typically 
up to a few tens per snapshot depending on the central concentration of the system. Given 
a set of blended stars, we only consider the brightest member of the ensemble as observable 
for the purpose of re-constructing the stellar mass function. As the luminosity of a main 
sequence star scales approximately with M 7 / 2 , this is an adequate approximation. With 
this procedure we define a catalog of observable main sequence stars where we include stars 
down to 0.2 M (for reference the observed mass fun ction of NCG 2298 h as greater than 



50% completeness at the center at this mass limit; see lPasquato et al.l 120091 ) . 



The stellar mass function in main sequence stars is modeled as a power law. Its index 
a is obtained by maximizing the likelihood L: 



( m a+l _ \ 
u d j+a lo ^ 

' i=l,N M SR 



(7) 



•MS.R 

where Nmsr is the number of resolved main-sequence stars and and mu are their mini- 
mum and maximum mass (here 0.2 and 0.8 M Q ). 

Based on the catalog of resolved main-sequence stars used above to retrieve a, we define 
the total observed main-sequence mass of the system M as: 

M= ^ m i- ( 8 ) 



3.2. 3-Dimensional, Mass-Based Structural parameters 



In a numerical simulation code such as NBODY6, the tidal radius {rt} s is defined based 
on Eq. [f5J We define as {M} s = is the total mass of the cluster (including dark remnants 
such as black holes, neutron stars and white dwarfs), that is the mass of all the particles 
within 2 {r t } s from the center of the cluster. 

Several different definitions have been used in numerical simulations for the core radius, 



so we direct the reader to lTrenti. Heggie fc Hutl (120071) for a com prehensive discussion. Here 
we adopt the standard definition from ICasertano &: Hutl (119851 ): 



(9) 



J2i=l,Ntot mi P i 

where is the distance from the center of the system and the sum is made over all the 
N tot particles of the system (including dark remnants). The density pi around each particle 
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is computed from the distance to the fifth nearest neighbor (ICasertano fc Hutlll985l ). In 
passing we note that this is not the N BODY6 definition of {r c } s , which is based on a p\ 
weight (see lTrenti. Heggie fc Hutll2007l ). The half-mass radius of a simulated system ({rh} s ) 
is defined as the tri-dimensional radius that contains half of the total mass of the system 
({M}J, thus including again any compact remnant that does not contribute the the surface 
brightness profile of the system. 



4. Results 

4.1. Mass Function Evolution 

The evolution of the mass function index a is shown in Fig. [T] as a function of the 
observed remaining mass of the system. As time progresses (and thus the system loses stars 
and mass), the mass function is preferentially depleted of low mass sta rs and becomes flatter 



(larg er a). Models in a stronger external tidal field lose star faster (ITrenti. Heggie fc Hut 



20071 ) and consequently have larger variation in a. However when looked at in terms of 
the remaining fractional mass, the evolution of a becomes independent of the tidal field 
strength and almost independent of the binary fraction considered in the runs. A larger 
binary fraction tends to moderately slow down the flattening of the mass function because 
low mass stars can be retained, and observed, if they have a dark remnant as a companion. 
All the simulations in Fig. [1] start with the same IMF below the turn-off (but with different 
retention fractions of neutron stars and black holes in the right panel). In addition, their 
initial conditions are such that the tidal field cut-off is self-consistent with the initial density 
profile of the system. To investigate if a universal evolution of the mass function index is 
present in a broader range of conditions, we consider in Fig. [2] additional runs that either 
start from a different IMF (Salpeter) or that are initially underfilling their Roche lobe. 
For the latter, we start from a Wo = 3 King model but we set {r t } s = 6.28, which is 
twice as large as the self-consistent {r t } s (Wq = 3) = 3.14. To compare the evolution of 
different mass functions, we consider the change in the power-law fit of the mass function: 
Act = a(t) — a(t = 0). From Fig. [2] it is immediately clear that the two simulations starting 
with a Salpeter IMF behave very similarly to the Miller&Scalo runs. The run initially 
starting with an underfilled Roche lobe does instead evolve marginally faster. Aa is in 
fact larger at a given fraction of mass loss, especially in the early stages of the evolution of 
the system. At later times the Roche lobe is eve ntually filled, consistent with the results 



of the simulations by iGieles fc Baumgardtl (120081 ). At that point the evolution of Aa is 



similar to our standard models. The close correlation between the slope of the stellar mass 
function and the fraction of the initial mass lost by a cluster was also found in previous 
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studies (jVesperini fc Heggid 119971 ; iBaumgardt fc Making 120031 ). Our simulations confirm 
this relation and show its validity for a broader range of different initial conditions. Recently, 



Kruijssenl (120091 ) investigated the evolution of the mass fun ction with an ana lytic model, that 



reproduces the universality of Aa for different IMF. The iKruijssenl ( 120091 ) conclusion that 
the retention fraction of black holes influences the evolution of the mass function (see Fig. 15 
and Sec. 5.2 in that paper) is howev er not seen in our simulations. This is likely because 
the analytic model of lKruijssenl (120091 ) does not take into account multiple encounters (three 
body or hi gher) of stellar black holes that contrib ute to thei r mutu al ejection on a relaxation 
timescale ( IMerritt et ajj 120041 ) . The models of IKruijssenl (120091 ) with different retention 
fractions start to show differences in the evolution of a when more than half of the mass 
of the system is lost. At that point only a few stellar mass BH would have survived in th e 
system, independently of their initial number (e.g. see lMerritt et al.l 12004 ; iGill et al.ll2008l ). 



The quasi-universal evolution of a as a function of MjMit = 0) is very interesting 
because, if th e initial mass f unction of Galactic globular clusters is universal, as for example 
suggested by iKroupal (120021 ) , then the slope of mass function is a direct tracer of the mass 
loss of the system, and thus of its dynamical age and of the tidal field it experienced during 
its life. This diagnostic is most effective to detect a large mass loss as \da/dM\ steepens when 
M decreases. Our conclusions on the quasi-universality of the a evolution can be applied 
to old globular clusters, with a turn-off mass below 1M Q . For younger systems both the 
turn-off mass and the ti dal disruption t imescale may introduce deviations from the universal 
behavior observed here (lKruijssenll2009l ) . Finally it is interesting to note that our conclu sions 
are consistent with Galactic globular clusters observations of a (IDe Marchi et al.l 120071 ): the 
clusters more massive than ~ 2 x 10 5 M© have essentially the same a, while lower mass 
clusters usually show a broader range of a values, likely indicating that mass loss has taken 
place. 



4.2. Structural parameters from King model fits 

The evolution of the observed structural parameters also tends to attain a universal 
configuration, as shown in Figs. [3] and H] for the core-to-half-light radius ratio and in Fig. [5] 
for the concentration c. The left panel of Fig. [3] shows the r c j rhp evolution of iV = 64k 
models with different tidal field configurations and initial density profiles as well as a binary 
fraction from 0% to 5%. Eventually all these models reach a common value for the observed 
fcj r hP ratio as derived from King profile fitting. However the observed r c j rhp ratio behaves 
differently from the theoretical and numerical expectation for {r c /rh} s - In fact, core-collapse 
for systems made only of single stars is detectable only for a very short time based on 
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^c/^hP- After core-collapse, runs with 10% primordial binaries, (left panel of Fig. H|) behave 
essentially as those with single stars only. We interpret this result as a consequence of 
mass segregation of dark remnants at the center of the cluster. The dark-remnants are on 
average heavier than the main sequence stars and thus they transfer kinetic energy to them. 
This results in a heating of the visible component of the cluster, which therefore does not 
exhibit a typical core-c ollapsed struc ture. An independent confirmation of this trend can be 
inferred from Fig. 5 in iHurleyl (120071) , wher e runs with and without primordial binaries also 
have a similar core (although IHurleyl 120071 runs do not reach past core-collapse). The core 
collapse for single-star systems is instead well discernible from the evolution of {r c /rh} s (see 
Fig. [6]). Until core-collapse, theoretical and observed King-fitted core radii evolve similarly, 
but immediately after the core bounce the two definitions separate sharply from each other. 
Post core-collapsed systems develop a shallow central cusp (/i ~ R~ u with v ~ 0.4 — 0.7) 
that is missed by the best fitting King model (see Fig. [7]), whose structure is constrained by 
the global surface brightness profile. In general, a King model becomes a poorer description 
of the system after core-collapse: in fact, the quality of our fits is degraded by about a factor 
two, as measured by the ratio of typical \ 2 before and after core-collapse. As a consequence, 
the fit results in the post core-collapsed phase also become somewhat dependent on the 
modeling of photometric errors. We assumed a constant fractional error on the surface 
brightness profile, which corresponds to a constant error in magnitude, except for the outer 
points, where we assigned a constant error in flux (to model the effect of the sky background). 
Assuming a constant error in flux for all the points in the surface brightness profile yields to 
no differences in the fit before core collapse. After the collapse, slightly more concentrated 
models are preferred, because in this case the points with the highest surface brightness have 
an increased weight relative to the outer points. 



The results in Fig. [3] have been obtained assuming a iMiller &: Scald (119791) I MF and a 
100% retention of neutron stars and stellar-mass black holes. A ISalpeterl (119551 ) IMF has 
a higher fraction of dark remnants per unit mass (because of it shallower slope, a = —2.3, 
compared to the MS IMF, a = —3.3, in the range [10 : 100]M o ), thus it sustains a larger 
core, especially in the earlier stages of the simulation (see right panel of Fig. PH). T his effect is 
primarily due to stellar mass b lack-holes (IMerritt et al.ll2004j ; iMackey et al.ll2008l ). so r c /ri t p 
returns in line with that of the IMiller fc Scald (Il979l ) simulations once the BHs segregate at 
the center of the cluster and kick each-other out of the system via three body interactions. 
Reducing the retention fraction of neutron stars and BHs to 30% leads to a decrease of the 
core radius down to r c /rh ~ 0.05, albeit with large fluctuations. The large fluctuations that 
appear after the core-collapse are also related to the degraded fit-quality, which increases 
the uncertainty on c (that is on W ). 



Interestingly, the main-sequence surface brightness profile has a significantly larger core 
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in the run with a central IMBH compared to all other runs once they r each their long-term 
equi librium configura tion (see right panel of Fig. H]). The prediction of iTrenti et al.l (120071 ) 
and iGill et al.l (120081 ) that collisionally rel axed systems with a ce ntral IMBH should exhibit 
a large r c /r hP is thus confirmed (see also iBaumgardt et al.ll2005l ). However systems with a 
large core radius might also be in the pre-core collapse stage. Unless the cluster relaxation 
time is short enough to suggest the cluster should have already undergone core collapse, no 
IMBH (or IMBH-like heating) is required to explain large values of r c /r/jp. A small value 
of r c /r hP , on the other hand, can be safely used to exclude likely candidates to harbor an 
IMBH. 

Fit of the surface brightness profile with a King model does not appear to provide a 
very solid determination of the intrinsic tidal radius of the system. The comparison between 
r t and {rt} s in Fig. [8] shows that, when the entire surface brightness profile is used for the fit, 
then r t > {r t } s . This is however not surprising, because unbound particle escaping from the 
system are removed from the calculation only when they have r > 2 {r t } s , thus the surface 
brightness profile at radii beyond {r t } s may still be positive. Similarly in actual star clusters, 
unbound stars remains initially in the proximity of the system, while escaping. The measure 
of r t becomes more uncertain, and dependent upon the details of the fit, if only a fraction of 
the surface brightness profile points is used (see blue dotted line in Fig. [8]). This is especially 
true after core-collapse (that is at M/M(0) < 0.65 in Fig. [8j This is primarily because we 
carried out the fit over the full-surface brightness profile by fixing the total luminosity and 
the half-light radius within 1% of the values in the profile. If only a fraction of the points is 
available, these two scales cannot be measured directly and must be left as free parameters 
of the fit. This leads to increased uncertainty, especially after core collapse, when the surface 
brightness profiles might present some departures from a King profile (for example see Fig. [9] 
at R ~ 0.15). When we do not use the outer points in the fit, the tidal radius tends to be 
moderately underestimated after core collapse compared to the determination using the full 
profile (Fig. IS]). This is consistent with wh at has been found in the determinatio n of the tidal 
radius for some Galactic globular clusters (IMcLaughling fc van der Marelll2005l ). A different 
concl usion may however be reached for those globular clusters that do not fill their Roche 
lobe IBaumgardt et al.ll20Q9h . 



4.3. Alternative core radius definitions 

Given the difficulty in reliably identifying core-collapsed systems from King model fit- 
ting, we explore here alternate definitions for the observed core radius. One possibility is 
to use the classical core radius r M , that is the radius at which the surface brightness falls to 
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half of its central value. A comparison between r c and r M is shown in Fig. [10] and suggests 
that this definition is better suited for identifying core-collapsed systems, remains in fact 
close to {r c } s (within ~ 15% accuracy), confirmin g for complex dynamic co nfigurations that 
include mass segregation the results published by ICasertano fc Hutl (119851 ) for King profiles 
with a constant mass-to-light ratio. However the definition of r M depends on the choice of 
bins, especially on the central one. Too small a bin size will lead to large Poisson fluctuations 
in the central surface brightness value, while an excessively large bin will overestimate the 
core size by lowering the central value of the surface brightness. A large central bin is not 
necessarily a choice, but may be imposed in actual obs ervations by a lim i ted an gular resolu- 
tion. To overcome some limitation related to binning, ICasertano Sz Hut Jl985h proposed to 
extend the density-based core radius estimator to the surface brightness profile f-i(R) and to 
define a surface-brightness density radius s^. In our notation this reads: 



E 



i=l,N MSj 



{i(Ri)Ri 



=1JV A 



(10) 



where the sum is carried out over all the detecte d main-sequence stars, s ,,. is shown for our 



simulations in Fig. [101 As expected from Fig. 3 in ICasertano fc Hutl (119851 ). s M deviates from 



r M and {r c } s for high-concentration King models, that is mainly in the post core-collapsed 
phase. Interestingly the value of r c obtained from King fitting for post core-collapsed systems 
is approximately bound between r M and s^. 



4.4. Fluctuations 



Fluctuations in the structural par ameters beyond those e xpected from Poisson noise 
are clearly apparent from our analysis. iHeggie &: Giersa (20091) reach t he same conclusion 
for their simulation of NGC 6397. Our run with a lMiller fc Scald (119791 ) mass function and 
low retention of dark remnants has particularly strong fluctuations, as can be seen for all 
different core radius estimators in Fig. [10J It is very interesting to note that the system 
appears to go through some breathing phases, where subtle changes are made to the surface 
brightness profiles, affecting primarily the King model fit. The origin of these fluctuations 
is likely dynamical, especially because we are considering only the main-sequence stars and 
thus we are relatively unaffected by low-number statistics. 
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4.5. Identification of core-collapsed systems 



Our results on the determination of core radii have profound implications for our un- 
derstanding of the dynamical state of observed globular clusters: most intrinsically core- 
collapsed globular clusters in our simulations appears to have a surface brightness profile in 
main sequence stars essentially identical to non-core collapsed systems when a King model is 
used to retrieve r c . Furthermore our simulations do not predict very small r c if the retention 
fraction of neutron stars and black holes is high. Alternative definitions of the core radius, 
such as r M , are more efficient at identifying core-collapsed clusters, but they may still suffer 
from large fluctuations and/or biases due to finite resolution at the center of the system. 
This picture clearly illustrates the potential problems in the observational identification of 
core collapsed clusters. In fact, the classification of an observed globular as core-collapsed 
or not might be a result of the combination of significant fluctuations in the total surface 
brightness profiles (enhanced with respect to our main-sequence on ly analysis because th e 
profile is defined primarily by a small number of giant stars; see also lHeggie fc Gierszl 120 09). 



comb ined possibly with the limited angular resolution of ground-ba sed data (ITrager et al. 



1995 ) used to construct the structural parameters presented in the Harris ( 19961 ) catalog. 
The core size in the post core- collapse phase depends mildly (oc l/loq(0.11N)) on the num- 
ber of particles in the system (jVesperini fc Chernoflill994j ; iHeggie. Trenti fc Hutll2006l ). thus 
it is expected that smaller core radii can be measured in systems with a larger number of 
particles. Further exploration of this issue is deferred to a follow-up paper that will include 
live stellar evolution to follow stars off the main sequence and thus to correctly model the 
total luminosity of the cluster. This will allo w us to carry out a more direct comparison with 
the ground-based data that define r c in the iHarrij (119961 ) catalog. 



4.6. Evolution in the (c; a) plane 



An interesting correlation has been recently reported by |Pe Marchi et al.l (120071 ) be- 
tween the slope a of the global stellar mass function of globular clusters and their central 
concentration c. Less concentrated clusters are more depleted in low mass stars than those 
with higher values of c. 

The evolut ion of a star clust e r in t he (c; a) plane is driven by collisional (two-body) pro- 
cesses, thus the |Pe Marchi et al.l (120071 ) observations test our understanding of the dynamics 
of globular clusters. As discussed in Sec. [Q, the observed correlation appears to be in conflict 
with the standard scenario globular cluster evolution which predicts that relaxation drives 
an increase in c as the cluster evolves toward core collapse wh ile at the sam e time preferen- 
tially depleting the system of low mass stars via evaporation (jSpitzerl 119871 ) . The evolution 
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of our simulated star clusters in the concentration-slope of the mas s function (c — a ) plane 
is shown in Fig. [TH where we also report the values observed by |Pe Marchi et all (20071 ) 
for those clusters that are well relaxed (t r h < 1 Gyr) according to the iHarrij (119961 ) cat- 



alog. In plotting the concen tration parameter c for the clusters s hown in Fig. [TT] we have 
used t he va lues obtained by iMcLaughling fc van der Marell (120051 ). which are based on the 



Kind (119661 ) models. Differences from the concentration values reported in |Pe Marchi et al. 



(200 7J) are however very smal l (Ac < 0.1, except fo r NGC 6712 which has c = 1.05 in 
McLaughling fc van der Marell fj2005h and c = 0.9 in be Marchi et all (boOTT ll While our 
simulations never reach a deep core-collapse (c > 2.5), their concentration tends to evolve 
toward c ~ 1.7, still larger than what observed for the systems with the most depleted stellar 
mass function. One possibility to explain the system with lower concentrations remains that 
of a very efficient source of heating at the center of the system, comparable with the heating 
generated by an IMBH (see cyan line in Fig. ITTT) . An IMBH itself appears to be excluded 
as the preferred explanation, because NGC 2298, where a central IMBH is rejected at high 
confidence level ( jPasquato et al.l 120091 ) . is one of the low c points in Fig. [TTJ Stellar mass 
black-holes are also very unlikely to be viable, because the heating source needs to remain 
efficient in the latest stages of the evolution of the system. Stellar mass BHs lose instead ef - 
ficiency at later times (see red line in the lower panels of Fig. [6} see also lMerritt et al.ll20 04) . 



It ha s been recently suggested that white dwarf kicks can sustain a large core (IFregeau et al. 



20091 ). However the effect of the kick has been shown only at the level of {r c /rh} s and thus 
its full impact on the observed concentration c is an open question. 

An alternative explanation for the discrepancy in the (c; a) plane evolution might rest 
in observational biases. While we carried out our analysis restricting to main-sequence stars, 
the tidal radius of observed globular clusters is mainly derived from ground-based data 
(ITrager et all If 9951 ). that thus are mostly dominated by giant-branch stars. Therefore mass 
segregation might lead to a systematic underestimate of r t from these observations compared 
to our analysis. 



Overall Fig. [TT] shows that the correlation between a and c observed by lDe Marchi et al 



(120071 ) is still in mild disagreement with the expectations from numerical simulations. Half 
of the 8 observed clusters are consistent with our models (including NGC 2298 which has 
(c = 1.31; a = 0.5) and could be explained by a core fluctuation). Two of them appear 
instead to be core-collapsed clusters. They have a conventional concentration c = 2.5, 
higher than our what we obtain from our simulations, as discussed in Section 14.21 (but they 
could be explained by assuming a low retention fraction of neutron stars and black holes). 
Two clusters have instead concentrations below our predictions, but here the consistency of 
simulations and observations might improve if data and models are analyzed with a fully 
self-consistent approach. 
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5. Conclusion and Discussion 

In this paper we discussed the evolution of the global mass function index a and of the 
structural parameters for multimass direct N-body models of star clusters in a galactic tidal 
field with and without primordial bi naries using a total number of particles up to N = 65536. 



We adopted a lMiller & Scald (119791 ) initial mass function evolved with an instantaneous step 
of stellar evolution to a turn-off mass of 0.8 M . The simulation results have been interpreted 
with great care in order to use definitions of structural quantities consistent with observations 
of Galactic star clusters. 

The initial density profile of the cluster or the tidal density field is not important in the 
long term, as the memory of the starting concentration is erased on a relaxation timescale, 
with a similar subsequent evolution for models starting from different initial conditions (see 
figs. CD to [5]). r c /rhp appears to be stable also over different IMFs, but its value depends on 
the retention fraction of neutron stars and stellar mass black holes. In particular a deep core 
collapse can be identified from the core radius as determined following the observational 
procedures, r c , only if few dark remnants are present in the system. Even a significant 
fraction (f=10%) of primordial binaries does not appear to influence much the observational 
determination of the core radius. This is likely related to the heating induced by mass 
segregation of dark remnants in the core of the system. 

These findings highlight that the core radius r c as obtained from King model fits over 
the complete radial extent of the surface brightness profile are unable to reliably identify a 
core-collapse cluster. Post-core collapsed systems are in fact relatively well fitted by medium 
concentration King models, although high residuals are present at the center of the sys- 
tem (see Fig. [7j). It is suggestive to note that these residuals appear similar to shallow 
cusps that have been obse rved in several galactic globular clusters with HST photometry 



( iNoyola fc Gebhardtll2006l ). In our simulations a shallow cusp (/i ~ R~ u with v ~ 0.4 — 0.7) 
in the surface brightness profile appears after core collapse (see Fig. [7]). The presence of such 
cusps may thus provide a hint that a system is past core-collapse, but other indicators are 
needed. One possibility if high angular resolution data are available, is to use as a measure 
of the core, that is the radius at which the surface brightness falls to half of its central value. 
r M does in fact track the evolution of the tri-dimensional core radius defined in simulations, 
i.e. {r c } s (see Fig. [TD]). In alternative, one could explore the use of the mass-to-light ratio if 
spectroscopic information are available. Shallow cusps at the center of medium concentration 
surface brightness profiles appears to be so common in the post core-collapsed evolution of a 
globular cluster that this indicator is not particularl y helpful in suggesting clusters likely t o 



harbor a central IMBH, contrary to past suggestions (IBaumgardt et al.ll2005l ; lMiocchill2007l ). 



Our investigation on the evolution of structural parameters "observed" from simulations 
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shows an improved agreement with the c — a correlation identified by lDe Marchi et al.l (120071 ) 
for Galactic globular clusters, but still the simulated clusters are too concentrated in their 
final stages of their evolution. This might be related to the definition of the tidal radius 
r t , and thus of th e concentration, from the > 15 years old, ground-based data used by 



Trager et al.l (119951 ). The light profile is in fact dominated by contributions from stars along 



the giant branch, more massive than average, and therefore more mass-segregated than the 
main-sequence stars we used in this analysis. In any case, from the observed c — a correlation 
we can expect, based on our analysis, that some of the systems with the most depleted stellar 
mass function can indeed be core-collapsed systems lurking in the sample of objects fitted 
by regular King profiles. This highlights the importance of carrying out a uniform survey 
to construct the structural parameters of Milky Way globular clusters, combining large area 
coverage with t he high angular-resolut ion data at the center, like those obtained for about 
50 globulars by ISarajedini et al.l (120071 ) with HST. 
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Fig. 1 . — Time evolution of the mass function index a for iMiller fc Scald (119791 ) runs as 
a function of the total observable mass remaining in the system M. Left panel: models 
with N to t = 65536 with different tidal fields and initial density profiles (King models with 
W = 3, 5, 7 — red, blue, black colors) and 0,2,5 % primordial binaries (solid, dotted, dashed 
lines). Right panel: Wq = 7 King models starting with N to t = 32768 particles and a binary 
fraction up to 10%. Black lines are a binary fraction 0-5%, the red line has f=10%. A model 
with no primordial binaries and a central IMBH is also shown in this panel (blue). The 
green line shows the model with 30% retention fraction of neutron stars and black holes. 
The evolution of the mass function slope is almost universal. 
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Fig. 2. — Time evolution of the mass function index a, shown as the difference from its 
t = value, Aa = a(t) — a(t = 0), for a series of runs with a Miller & Scalo and a Salpeter 
mass function. The black line refers to a Miller & Scalo run with 100% retention fraction 
of dark remnants and no binaries starting from Wo = 7. The blue line has the same IMF 
but starts a Wq = 3 King model and {r t } s = 6.28 (this model has its Roche lobe underfilled 
by a factor 2). The red lines are runs with a Salpeter IMF, starting from Wq = 7 and no 
binaries: (solid: 100% retention fraction, dashed: 30% retention fraction). All simulations 
have N tot = 32768. 
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Fig. 3. — Evolution of the observed core to half light radius for our N = 64k runs, color 
codes like in the left panel of Fig. [U Left panel shows r c /rhp as a function of M, right panel 
as a function of t in units of the initial half-light relaxation time. 
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Fig. 4. — Evolution of r c /r^p like in Fig. [3] but for models with N = 32k. Upper set: 
Miller & Scalol Jl979h IMF and different binary fractions (0% black, 1% blue, 3% yellow, 



5% magenta, 10% r e d). L ower set: models with no binaries but d ifferent IMFs. Black 
line: iMiller fc Scalol (119791 ) with full retention of remnants, red line ISalpeterl (119551 ) with 
full retention of remnants. Cyan and magenta lines are the corresponding IMFs with 30% 
retention fraction. Blue is our run with a central IMBH. Only a central IMBH is able to 
sustain a large r c in the long term, but a significant number of stellar-mass BHs can also 
give a transient high value of r c . Small r c values are instead possible with a low retention 
fraction of dark remnants. 
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M/M(0) t/t rh (0) 

Fig. 5. — Evolution of the observed concentration c for our N = 64k runs, color codes like 
in Fig. [U Left panel shows c as a function of M, right panel as a function of t in units of 
the initial half-light relaxation time. 
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Fig. 6. — Observed core to half-light radius ratio (r c /rhp solid black line) and 3D-density 
based core to half-mass radius ratio ({r c /r h } s dashed blue line) for a W = 7 model without 
primordial binaries (N=64k, upper left and N=32k, upper right) as well as with 3% binaries 
(N=32k, bottom left) and with 10% binaries (N=32k, bottom right). The two definitions 
differ significantly after core collapse, especially for models with a low fraction of primordial 
binaries. Core collapse signatures are clearly present in {r c /rh} s for models without binaries 
but not in r c /rhp- 
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Fig. 7. — Surface brightness profiles (arbitrary units for /i, NBODY units for R) for eight 
snapshots of our iV = QAk, Wq = 7 simulation with no binaries. The snapshots are taken 
at times t = 470, 2020, 2660, 2940, 3000, 4000, 6000, 8000 in code units (equivalent to 
~ 1.0, 4.4, 5.8, 6.4, 6.5, 8.7, 13.1, 17.4 t r h(0)) and represents a typical pre-core-collapse 
profile, a profile on the way to core-collapse, the profile at the peak of core-collapse and five 
post-core-collapsed profiles. Superimposed to the data (black crosses) the best fit King model 
is shown. The five models have W = 7, 8.1, 9.4, 7.7, 8.0, 7.9, 7.7, 7.7. The quality of 
the fit in the post-core-collapse profile is about two times worse than in the other snapshots. 
In the post-core collapse profiles (last five panels) a shallow central cusp (/i ~ R~ u with 
v ~ 0.4 — 0.7) is present (blue dashed line). 
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Fig. 8. — Tidal radius as measured observationally using the full surface brightness profile 
(r t — solid black line) or only the inner part containing 75% of the total light (dotted blue 
line) in our N=64k Wq = 5 simulation with no primordial binaries. If information on the 
outer parts of the system is missing, then the tidal radius determination becomes dependent 
upon the details of the fitting method after core-collapse (M/M(0) < 0.65). The theoretical 
tidal radius is also shown {{r t } s — red dashed line). At later times r t j {r t } s ~ 2 when using 
the full surface brightness profile for the King model fit. 
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Fig. 9. — Surface brightness profile (blue crosses) at t=3420 (in the post-core collapse phase) 
for iV = 64k simulation with no primordial binaries starting from Wq = 5. The standard 
King model fit over the full radial range of the profile, and with total luminosity and half- 
light radius fixed to the measured quantities, is shown as black solid line (Wq = 7.5). A King 
profile fit using only the inner 75% of surface brightness profile points and with free total 
luminosity and half-light radius is shown as red solid line (Wo = 7.3). The dotted vertical 
lines (black and red) show the outer boundary of the region used for fitting the King models 
in the first and second case. 
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Fig. 10. — Comparison of different core radii definitions for two of our simulations without 
primordial binaries (left panel, N=64k, Wo = 5; right panel N=32k, Wq = 7, retention 
fraction 30%). Post core-collapse oscillations in r c are particularly strong in the right panel. 
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Fig. 11. — Mass function index a versus the central concentration for the iV = QAk models of 
fig. [Hand for th e N = 32k mode l with a central IMBH (cyan solid line). The green squares 
are points from |Pe Marchi et al.l (120071 ) associated to globular clusters with t r h < 1 Gyr. 
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Table 1: Summary of N-body simulations 



AT 

iV 


J 


TTV /TT~1 

Mr 


T XT 


65536 


(J. (J (J 


A /TO 


6 


65536 


0.02 


MS 


3 


65536 


0.055 


MS 


3 


65536 


0.00 


MS 


5 


65536 


0.02 


MS 


5 


65536 


0.00 


MS 


7 


32768 


0.00 


MS 


7 


32768 


0.01 


MS 


7 


32768 


0.03 


MS 


7 


32768 


0.05 


MS 


7 


32768 


0.10 


MS 


7 


32768 


0.00 


Sa 


7 


32768 


0.00 


MS 


7 


32768 


0.00 


Sa 


7 


32769 


0.00 


MS 


7 


32769 


0.00 


MS 


3 



Other 



30% NS/BH retention 
30% NS/BH retention 
miMBH = 0.01 
{r t } s = 6.28 [underfilled 2] 



Note. — N-body simulations of star clusters in a tidal field with self consistent King model initial con- 
ditions. First column: number of particles N; second column: binary fraction / = Nb/N; third column: 
initial mass function used (Sa: power law, MS: Miller & Scalo); fourth column: initial concentration of the 
density profile (King index Wo)- Fifth column reports other relevant information for the simulation. The 
last simulation starts from compact initial conditions where the Roche lobe is underfilled by a factor 2. 



